arXiv: 1506.08281 v3 [cond-mat.mes-hall] 20 Nov 2015 


Ah initio analysis of the topological phase diagram of the Haldane model 


Julen Ibanez-Azpiroz,^ Asier Eiguren,^’^ Aitor Bergara,^’^’^ Giulio Pettini,^ and Michele Modugno®’^ 

^Peter Griinberg Institute and Institute for Advanced Simulation, Forschungszentrum Jillich & JARA, D-52425 Jiilich, Germany 
^Depto. de Fisica de la Materia Gondensada, Universidad del Pais Vasco, UPV/EHU, 48080 Bilbao, Spain 
^ Donostia International Physics Genter (DIPC), 20018 Donostia, Spain 
'^Centro de Fisica de Materiales GFM, Gentro Mixto GSIG-UPV/EHU, 20018 Donostia, Spain 
^Dipartimento di Fisica e Astronomia, Universitd di Firenze, and INFN, 50019 Sesto Fiorentino, Italy 
^Depto. de Fisica Teorica e Hist, de la Ciencia, Universidad del Pais Vasco UPV/EHU, 48080 Bilbao, Spain 
^IKERBASQUE, Basque Foundation for Science, 48011 Bilbao, Spain 
(Dated: November 23, 2015) 

We present an ab initio analysis of a continuous Hamiltonian that maps into the celebrated 
Haldane model. The tunnelling coefficients of the tight-binding model are computed by means of 
two independent methods - one based on the maximally localized Wannier functions, the other 
through analytic expressions in terms of gauge-invariant properties of the spectrum - that provide 
a remarkable agreement and allow to accurately reproduce the exact spectrum of the continuous 
Hamiltonian. By combining these results with the numerical calculation of the Chern number, we 
are able to draw the phase diagram in terms of the physical parameters of the microscopic model. 

Remarkably, we find that only a small fraction of the original phase diagram of the Haldane model 
can be accessed, and that the topological insulator phase is suppressed in the deep tight-binding 
regime. 

PACS numbers: 67.85.-d,73.43.-f 


I. INTRODUCTION 

The Haldane model [T] is a celebrated lattice model de¬ 
scribing a Chern insulator [5], characterized by the pres¬ 
ence of quantum Hall effect (QHE) in the absence of 
a macroscopic magnetic field. Conceptually, the Haldane 
model stands at the heart of the tremendous advances 
in the field of topological condensed matter physics, as 
the mechanism for a non-trivial band topology presented 
by Haldane is realized in actual materials via the intrin¬ 
sic spin-orbit interaction of topological insulators nig. 
These concepts are also relevant for the physics of ultra¬ 
cold atoms in optical lattices, as these systems represent 
a powerful platform for simulating solid-state physics [6] . 
Mott insulators [HIE], bosonic superfluids [3] or graphene¬ 
like honeycomb lattices [TUHTB] are among the many sys¬ 
tems that have been emulated by this technique. In¬ 
terestingly, an effective experimental realization of the 
Haldane model has been recently reported in Ref. HU. 

In his original work, Haldane constructed a discrete 
tight-binding model for a non-centrosymmetric honey¬ 
comb lattice in the presence of a vector potential A{r), 
with vanishing total flux through the unit cell. The key 
feature of the model is that, even in absence of a macro¬ 
scopic magnetic field, the time-reversal symmetry is bro¬ 
ken due to the presence of the gauge held A{r). This, 
in turn, implies that the next-to-nearest neighbour tun¬ 
nelling coefficient ti becomes a complex number. Hal¬ 
dane showed that the properties of the system depend on 
the interplay between the phase acquired by ti and the 
effect of parity breaking, affecting the topological phase 
diagram of the model [1]. 

Considering the above, the knowledge of the depen¬ 
dence of the phase acquired by ti on the applied vector 


potential held becomes a crucial element for drawing the 
topological phase diagram. For this purpose, it is com¬ 
mon practice [DUB] to make use of the so-called Peierls 
substitution, whereby the effect of A{r) is effectively in¬ 
cluded by the replacement ti ti exp {i{e/K) J A{r)dr) 
[El- However, in a recent work [20] we showed that the 
Peierls substitution is actually wrong whenever the vec¬ 
tor held A{r) has the same periodicity of the underlying 
lattice, as it is the case of the Haldane model by construc¬ 
tion. In that work, we analyzed the parity invariant case 
by presenting two independent approaches for calculat¬ 
ing the tight-binding parameters of the model: one based 
on the maximally localized Wannier functions (MLWFs), 
the other on a closed set of analytical expressions in terms 
of the energy spectrum at selected high symmetry points 
in the Brillouin zone (BZ). 

In the present work, we extend the previous analysis 
to the general case in which both inversion and time- 
reversal symmetry can be broken. We show that the two 
approaches considered provide a remarkable agreement 
even in the presence of parity breaking, allowing for a 
precise determination of the tight-binding parameters of 
the model. By combining these results with the numer¬ 
ical calculation of the Chern number, we are able to re¬ 
draw the topological phase diagram of the Haldane model 
in terms of the physical parameters of the microscopic 
model. Interestingly, we find that only a small fraction 
of the original phase diagram can be accessed, and that 
the topological insulator phase shrinks dramatically as 
the system becomes more and more tight-binding. In 
addition, we find that the gap closing at the topological 
phase transition does not take place exactly at one of 
high symmetry points of the BZ, but in a close-by point. 
The reason is that the complex tunneling between ho- 
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mologous sites are no longer degenerate in the presence 
of parity breaking, contrarily to what it is assumed in the 
Haldane model. 

The paper is organized as follows. In section [TT| we in¬ 
troduce the microscopic continuous Hamiltonian used in 
this work and review the formal steps needed to derive 
the corresponding tight-binding model. Some general 
properties of the Haldane model are also recalled. Then, 
in section |HI| we present the two approaches employed 
for calculating the tight-binding parameters, and discuss 
how the breaking of time-reversal and/or parity affects 
their behavior. In section [TV] we analyze the topological 
phase diagram, both in terms of the parameters of the 
tight-binding model and of the physical ones. Concluding 
remarks are drawn in section |V] Finally, in the appen¬ 
dices we present an analysis of the spread functional of 
the MLWFs (Appendix]^ and additional remarks on the 
numerics (Appendix [b| . 


II. SETUP OF THE HALDANE MODEL 


(shown in Fig. is generated by the direct lattice vectors 
ai ,2 = (27r/3fcL)(e2,=F'\/3ey). We also define the lattice 
vector = ai + a 2 , that will be useful later on. 

We now turn to the vector potential contained in the 
Hamiltonian ([I). As already mentioned, we employ the 
form proposed by Shao et al [T5], namely 


A{r) =ahki^ 


(sin ((62 


bi) ■ r)+ 


^^(-l)*sin(bi • r))e^ 


■ r)ey 


2=1 


(3) 


with V • A{r) = 0 (Coulomb gauge). The flux of the cor¬ 
responding magnetic held B — V x A across the unit cell 
is null [18] , as required for the Haldane model. In the fol¬ 
lowing analysis, the only variable parameter entering the 
expression for the vector potential A{r) is its amplitude 
a. Notice that for a = 0 the system is symmetric under 
time-reversal, whereas this is not the case for a ^ 0 . 


In this section we present a systematic derivation of the 
Haldane model starting from the continuous Hamiltonian 
proposed by Shao et al. HHI in the context of cold atoms 
trapped in optical lattices (see also HO]). The method 
discussed here is general and suited to map a generic con¬ 
tinuous Hamiltonian to its corresponding tight-binding 
model HUHO]. 

A. The continuous Hamiltonian 


B. The tight-binding model 

The continuous Hamiltonian Q can be mapped onto 
the tight-binding Haldane modeljll [18] by following the 
general procedure discussed in HOHOl]- The starting 
point is the (single particle) many-body Hamiltonian de¬ 
fined by 

-Ho = j dr (4) 


The general form of a continuous Hamiltonian in the 
presence of a scalar lattice potential VL{r) and a vector 
potential A{r) is 

Ho = ^[P-A{r)f+ VL{r), ( 1 ) 


with r = (x, y) in case of a two dimensional system, as 
we shall consider here. The potential VL{r) is chosen in 
order to generate a periodic structure with two minima 
per unit cell, forming a honeycomb lattice HnillSl: 


VL{r) = 2 sEr^ cos [( 6 i - 62 ) • r] 

+ cos {bi-r- + cos (62 • r) |. 


( 2 ) 


Above, Eji = Ti^k\/2m is the recoil energy, de¬ 
notes the laser wavevector, s is a dimensionless param¬ 
eter representing the strength of the potential in units 
of Er, bi ^2 = {V^kL/2.){ex'^^/iey) are the basis vec¬ 
tors in the reciprocal k space, and xa is a parameter 
related to the breaking of the parity symmetry. In par¬ 
ticular, XA = 0 corresponds to the inversion symmetric 
case, where the two minima in the unit cell are degen¬ 
erate. On the other hand, for xa ^ 0 parity is broken 
and the minima are no longer degenerate. The unit cell 


with a field operator for bosonic or fermionic particles. 
Then, when the wells of the lattice potential are deep 
enough, the held operator can be conveniently expanded 
in terms of a set of functions Wj^{r) localized at each 
minimum: 


= (5) 

jb' 

Above, ^ is a band index and (%i/) is the creation 
(destruction) operator of a single particle in the jth cell, 
satisfying the usual commutation (or anti commutation) 
rules following from those of the held ip. 

As already mentioned in the introduction, we shall 
use the maximally localized Wannier functions (MLWFs) 
for composite energy bands [53] as basis functions. The 
MLWFs are dehned as linear combinations of the Bloch 
eigenstates ip,y'k{'r'), 

If . ^ 

= -7^ / dk ^ U^ui{k)ip^.k{r), ( 6 ) 

where Hg represents the volume of the hrst BZ, and 
U G U{N) is a gauge transformation that is obtained 
from the minimisation of the Marzari-Vanderbilt spread 
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functional, — {'>')V\ US]- We remark that 

the presence of the vector potential may significantly af¬ 
fect both the Bloch eigenfunctions and the uni¬ 

tary matrices Uy^i{k) [5D]. A thorough analysis of the 
MLWFs is given in section [III A| 

In the following, we shall consider the contribution 
of the first two Bloch bands only, namely v,v' = 1,2. 
This is sufficient for constructing the lowest lying ML¬ 
WFs localized at the two lattice sites A and B in¬ 
side the unit cell. Then, by considering the follow¬ 
ing transformation from coordinate to reciprocal space, 
buk = CLjy, the reduced tight- 

binding Hamiltonian can be written as 

^ f dkK,,{k)bl,^K,k, (7) 
where the 2x2 matrix (fe) is given by 

huv'{k) = {wou\Ho\wjvi) 

3 

= f U*^{q)U,,>n{q)en{q), 

j -ISts 

( 8 ) 


with Rj = {jiSii -k j 2 a 2 ji,j 2 = 0,±1,±2 ...} and j 
labels the unit cell. We remark that the eigenvalues 
of h^^i{k) coincide with the exact Bloch energies e„(fe) 
(u = 1,2) if the full expansion of neighbouring coeffi¬ 
cients Rj is retained. When the system is in the tight- 
binding regime (s ^ 5) |21j . it is convenient to truncate 
the series by retaining only a finite number of matrix 
elements {wo^\Ho\wj,j'), with the eigenvalues of h^^'{k) 
still being a good approximation of the exact energies. 
We note that since the functions Wj^{r) are in general 
complex (see section III A), we may expect the matrix 
elements to be complex as well. 

By truncating the tight-binding expansion in Eq. Q 
to the next-to-nearest neighbours (see Fig. [^, the matrix 
huv'{k) can be written as the sum of three terms: 


K.'{k) = [h^°i{k) ++/iLlUfc)- (9) 

The first term corresponds to the on-site energies 

h^!^}{k) = Ey = {wq^\Hq\wo^), (10) 

which are real quantities by definition. The next term, 
contains only off-diagonal elements corresponding 
to the hopping between the three nearest-neighbour sites. 
Although the basis functions Wju{r) are complex, these 
three tunnelling amplitudes can be chosen to be real by 
means of a suitable global gauge fixing, as they are all 
equal thanks to the symmetries of the system (see Fig. 
[^. Taking this into consideration, and defining 

( 11 ) 



FIG. 1. (color online) Bravais lattice associated to the hon¬ 
eycomb potential in Eq. (|^. Black and white circles refer to 
minima of type A and B, respectively. The elementary cell 
is highlighted in grey. The various tunnelling coefficients are 
indicated for the site of type A in the central cell. The sys¬ 
tem is invariant under discrete translations generated by the 
Bravais vectors ai /2 and under rotations oi 6 = 27r/3 radians 
around any vertex of the lattice. The rotational symmetry 
implies that next-to-nearest tunnelling amplitudes along 
the same direction are conjugate pairs (solid and dashed lines 
in red); from the latter follows the equivalence of the hopping 
amplitudes separated by 27r/3 radians. When sites A and B 
are degenerate, the system is also invariant under rotations 
by TT radians around the centre of any elementary cell. 


we can write 


h^llik) = to{l + = toZoik) = z{k). 

( 12 ) 

Its conjugate counterpart is given by hW (fc) = z*{k). 
(‘2') 

Finally, the term hlJ (k) corresponds to the six next-to- 
nearest neighbours hopping between homologous sites. 
By taking into account all the symmetries of the full 
Hamiltonian of Eq. Q, these tunnelling coefficients can 
be compactly written as 

= {woi,\Ho\w±aj^) = , j = 1,2,3. (13) 


The above equation explicitly shows that contains 
two distinct complex phases for each site type {u = 
A,B). Then, using Eq. (13) and after some algebra, we 
can write 


hl^Xk) =\t iy||2cos [k as + + 2 ^ cos (fc-a^ — ip^) | 


i=l,2 


= \ti^\F^{k) = f^[k). 


(14) 


The above expressions allow to cast the matrix (fe) 
in the following compact form. 


feiyiy'(fe) — 


to = (woa|77o|wo_b) 


e^(fe) z{k) 
z*{k) tB{k) 


(15) 
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where we have defined 

e„{k) = E^ + U{k). 


(16) 


By expanding the Hamiltonian on the basis of the 2x2 
identity matrix, /, and of the Pauli matrices, tJi, Eq. (15) 
can be rewritten as [15] 


h{k) = ho{k)I + ^ hi{k)(Ti, 


(17) 


i=l 


where the coefficients hi{k) are given by the following 
expressions: 

eA{k) + eB{k.) fA{k) + fB{k) 


hoik) = 

= /+(fe), 

3 

/ii(fe) = Re[z(fe)] = to ^ cos (fc • Si ), 

i=l 

3 

/i 2 (fc) = -Im[z(fc)] = to ^ sin (fc • s*), 


i=l 


i. eAik)-eB{k) 

hoik) = - - -= e - 

= e + /_(fc). 


fAik)-fBik) 


(18) 

(19) 

( 20 ) 
( 21 ) 


with the vectors Si being those joining the three nearest 
neighbours of type AiB) to a given site of type BiA) [TB]. 
In the last expression we have also fixed, without loss of 
generality, Ea = e and Eb = — £■ 

Then, the tight-binding energies are readily found as 

e±(fc) = hoik)±\hik)\ = /+(fc)±\/[e + /-(fc)]^ + k(fc)P, 

( 22 ) 

where h = (fci, h 2 ,ho). 


1. The Haldane model and the Peierls substitution. 


At this point, further approximations are required in 
order to recover the original model proposed by Haldane 
[I], namely |ti^| = \tiB\ = \ti\ and tpA = -pB = T- We 
shall refer to this configuration as the “simplified param¬ 
eter setup” (SPS). We note that the corresponding model 
contains only four parameters, namely e, toi |^i| and ip. In 
section [III C[ we shall provide numerical evidence show¬ 
ing that in the tight-binding regime the difference be¬ 
tween |ti^| and [tisl is negligible, thus justifying the 
SPS. The second condition is not strictly verified (again, 
see section HIC), but one can consider a sort of effec¬ 


tive model by defining a single phase, ip = ip a — 'Pb)/ 2. 
Therefore, in the SPS the terms ho and h^ of Eqs. (18) 
and (211 simplify to 


ho = 2\ti \ cos sin (fc • a^) , 

i=l 

3 

fc 3 = e — 2|ti| siny:^^ sin (fc • Oi). (23) 


The above equations, together with Eqs. (19) and 


(20), correspond to the definition of the original Haldane 
model, see [B US¬ 
As already mentioned in the introduction, the original 
model proposed by Haldane is constructed by means of 
the so-called Peierls substitution BilH!- This consists 
in assuming that the complex phase of the tunnelling 
coefficient tij is given by the integral of the vector po¬ 
tential along the straight path joining sites i and j, i.e. 


ti 


t 


exp{iie/h) // Adr}. 


In the present case, the 


Peierls prediction for the complex phase would be 


27r 


(24) 


This value will be used later on, in Secs. HI and IV 


for comparison with the results of the two approaches 
discussed in this paper. Here, we can anticipate - as 
thoroughly discussed in [50] - that the above result is 
definitely incorrect, owing to the fact that the Peierls 
substitution is justified only when the vector field A(r) 
is slowly varying on the scale of the lattice [53|. In fact, 
this condition is explicitly violated in the Haldane model, 
where A(r) has the same periodicity of the lattice. 


C. General features of the Haldane model 

Before proceeding, let us recall some general features 
of the Haldane model [Tj, corresponding to the SPS ap¬ 
proximation. This model is characterized by the presence 
of Dirac points located at the vertices kjj of the first 
BZ, where the dispersion law takes the relativistic form 
e(fc) = \/rrBA + fc^. They can be divided into two in¬ 
equivalent classes, corresponding e.g. to fc^ = ±(l,0)fcL 
(often referred in the literature as K and K'; in the follow¬ 
ing, we will always refer to these two as inequivalent Dirac 
points, for simplicity). In the presence of time-reversal 
and inversion symmetry (namely, for a = 0, xa = 0)) the 
two energy bands become degenerate at the Dirac points, 
whose position is given by the solution of ^(fcu) = 0 or 
hi(fcz)) = h 2 (fc_D) = 0 (this holds at any order of the 
tight-binding expansion). In fact, in this case both e and 
p vanish, yielding hoik) = 0 , /a = /s, and /_ = 0, so 
that the energies e±(fc_D) are degenerate. 

When both time-reversal and inversion symmetry are 
broken (a ^ 0 ,Xa 0)) two inequivalent energy gaps 
form at the Dirac points: 

= e+(fc±) - e_(fc±) = 2\hoikU = 2y[e + /_(fc±)]2. 

(25) 

The closure of one of them indicates a topological phase 
transition, where 


6 += 2 


:±3y3|t 


1 1 Sin (f 


= 0 . 


(26) 


This equation identifies the well known boundary be¬ 
tween the normal and topological insulator phases with 
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Chern numbers (7 = 0 and (7 = ±1 in the Haldane model 

m- 

We remark that in the general model the gap closing 
does not take place exactly at but in a close-by point. 
This is due to the fact that, when breaking of parity is 
included self consistently, the tunnelling parameters tiA 
and tiB are no longer degenerate, contrarily to what it is 
assumed in the Haldane model (cf. the SPS approxima¬ 
tion). 


III. CALCULATION OF THE TIGHT-BINDING 
PARAMETERS 

In this section we discuss two independent methods 
for calculating the tight-binding parameters for arbitrary 
values of the physical parameters s,a and XA- The first 
method is based on the ab initio calculation of the max¬ 
imally localized Wannier functions (MLWFs) [231 US], 
which we already employed in I20H22] in different lat¬ 
tice geometries. This approach gives direct access to the 
whole set of parameters e, tp: and ips- The 

second approach relies instead on analytical expressions 
in terms of the energy spectrum, as discussed in [20j . The 
latter is here extended to the case of parity breaking. 

We remark that the approach based on the MLWFs 
corresponds to the ab initio definition of the parameters, 
and is therefore model-independent. Instead, the second 
method depends on the specific form of the tight-binding 
Hamiltonian. However, it does not require the calcula¬ 
tion of any set of Wannier functions since only the spec¬ 
trum of the continuous Hamiltonian is needed. Notably, 
the two methods present a remarkable agreement in the 
whole range of parameters considered here. 


A. Maximally localized Wannier functions 

The MLWFs, which have been defined in Eq. ®> , rep¬ 
resent a powerful tool that is largely employed in con¬ 
densed matter physics |2S|. By construction, the ML¬ 
WFs are the basis functions with the maximal degree 
of localization in coordinate space, allowing to construct 
tight-binding models that accurately reproduce the prop¬ 
erties of the continuous Hamiltonian, with a minimal set 
of tunnelling coefficients. In addition, the MLWFs per¬ 
mit a very fine sampling of the reciprocal space thanks to 
the so-called Wannier interpolation technique |2S]. This 
point is very important for our purposes in this work, as 
the determination of the Chern number requires a high 
density of points in fc-space [Ml El] ■ 

The MLWFs are computed by means of the standard 
implementation of the WANNIER90 package [Mj (see 
also Appendix]^. The resulting functions are complex¬ 
valued when a ^ 0. This feature is in agreement with 
the analysis of [29], where it was shown that, in general, 
MLWFs cannot be constructed as real functions when 
the time-reversal symmetry is broken (see also [30] 1. In 


(a) log|Re[wJr)]f 



o’ ‘ ‘ ‘ ‘ ^.2 ‘ ‘ ‘ ‘ 3 


x(units of a) 

FIG. 2. (Color online) Density plot (in logarithmic scale) 
of the square of the real (a) and imaginary (b) parts of the 
MLWF for sublattice A, for s = 5, a = 0.1 and xa = 0. The 
solid and dashed lines denote the unit cell and the honeycomb 
lattice of the scalar potential, respectively. In the latter, the 
corners of the hexagons mark the minima of the scalar lattice 
potential labelled either as A or R. 


the context of the Haldane model, the imaginary part of 
the MLWFs plays an essential role, since it determines 
the complex phase acquired by the next-to-nearest tun¬ 
nelling coefficient. In turn, the complex phase directly 
affects physically meaningful quantities, such as the spec¬ 
trum or the topological phase diagram [T]. 

In Fig. [^ we illustrate an example of the real-space 
structure (note the logarithmic scale) of the real and 
imaginary parts of a MLWF for sublattice v = located 
at the origin j = 0, for s = 5, a = 0.1 and xa — 0. The 
structure of the real part is very similar to the one in the 
pure honeycomb lattice mi, namely it is highly localized 
around the origin, with appreciable contribution around 
the neighbouring lattice sites. In average, the imaginary 
part is two-three orders of magnitude smaller than the 
real part. It is particularly interesting to observe that the 
imaginary part is null at the interstitial region between 
nearest neighbours, while it becomes maximum along the 
path joining next-to-nearest neighbours. These proper¬ 
ties hold in the whole range of parameters considered in 
this work. 
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An analysis of the spread functional of the MLWFs 
as a function of the vector potential amplitude has been 
included in Appendix 


B. Analytical expressions from the spectrum 


In this section we derive a closed set of analytical ex¬ 
pressions in terms of the energy spectrum at selected high 
symmetry points in the BZ. This is done in the framework 
of the SPS discussed in section IIB ![ corresponding to 
the standard formulation of the Haldane model fflUH]. 
As we shall see below, the approximations of the SPS 
are well justified in the tight-binding regime. The model 
is therefore given in terms of four parameters, namely 
e,(p,to and |ti|. We remind that e measures the differ¬ 
ence between the on-site energies and Eg, and it is 
therefore associated to the breaking of parity, whereas 
the breaking of the time-reversal symmetry corresponds 
to if different from zero. It is also worth recalling that 
the parameters of the underlying continuous Hamiltonian 
that control the breaking of parity and time-reversal sym¬ 
metry are XA and a, respectively. In particular, XA = 0 
gives e = 0 whereas a = 0 implies (p = 0. 

We begin by noting the following relations at fc = 0; 


the following set of formulas: 

^+ + <5- 


1 


^0 “ r\/(^++ +'^+) “ 


2 (5+ -I-15_) 




(p = tg 


-1 


■\/3 (5+ — S- 
2 A+ - 


(37) 

(38) 

(39) 

(40) 


Similarly, in the region with e < 3-\/3|ti| sin:/? (corre¬ 
sponding to the topological insulator phase), we find the 
following expressions: 




1^1 1 = Yg\/(^+-^-) +^(^+ + ‘5-)^. 


(fi = tg- 


73 5+-h(5- 

^ A+-At 


(41) 

(42) 

(43) 

(44) 


/+(0) = 6|ti|cosv?, 
/-(O) = 0, 

170)1 =3|to|. 


(27) 

(28) 
(29) 


The solutions in a generic case with e < 0 or (/? < 0 can be 
obtained from symmetry considerations, by exchanging 
the role of the two basis points A, B and/or of the two 
inequivalent Dirac points 


Similarly, at the Dirac points k 


D 


we have 


C. Numerical results 


f+iko) =-^\ti\cosip, (30) 

/_(fe^) = ±373|ti| sinv?, (31) 

Next, let us define the bandwidths 

A±=+[e+(0)-e+(fc±)], (33) 

A^ = -[e_(0)-e_(fc±)]. (34) 


Recalling the expression for the gap at the Dirac points 
in Eq. (26), one can easily derive the following relations: 



A+ -H At -h h+ 
2 


A^_ -|- A_ -|- S— 
2 


18|ti I cos (f = A/i — At = A_|_ — A_. 


,(35) 

(36) 


Due to the symmetries of the system, we can consider 
e > 0 and (p > 0 without loss of generality. Focusing first 
on the region with e > 3-\/3|ti| sin:/? (corresponding to 
the normal insulator phase), after some algebra one finds 


In this section we present a comparison of the two 
methods described in Secs. IIII Al and IHIBI for the cal¬ 
culation of the tight-binding parameters. In addition, we 
also analyze the accuracy of the assumptions of the SPS 
(Sec. |!h^ based on the tunneling coefficients extracted 
from the MLWFs. 

Let us begin by analyzing Fig. where we compare 
the tunnelling coefficients calculated from the MLWFs 
with those calculated from the analytical formulas of 
Eqs. (37)-(40), valid for the normal insulator regime, and 
Eqs. (411-(44), valid for the topological insulator regime, 
which is depicted by the grey shaded area in the figure. 
Results are shown as a function of a for fixed values s = 5 
and XA = 0.001, since essential features are unaffected by 
s and XA- In the case of the MLWFs, we have plotted the 
averages (p = (ipA - and |ti| = (Iti^l + |tiB|)/2 in 

order to allow comparison with the analytical formulas, 
which have been derived in the context of the SPS (see 
section 


HIB) 


Overall, Fig. shows a very good agreement between 
the two methods for all the tunnelling coefficients, in 
all regimes. Furthermore, it is interesting to note that 
the two different solutions represented by the set of Eqs. 
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a 


FIG. 3. (color online) Comparison of the four tight-binding 
coefficients as calculated from the MLWFs (solid red lines) 
and the analytical formulas of Eqs. ( |37[ )-(|40| (blue triangles) 
and Eqs. (41l-(44l (green circles). Results are shown as a 
function of tr, keeping fixed values s = 5 and xa = 0.001. 
The grey area in the figures denotes the region where the 
system behaves as a topological insulator with C 7 ^ 0 (see 
text and section |IVI) . 



0.0 0.5 1.0 1.5 2.0 


a 


FIG. 4. (color online) Relative deviations from the average 
values of (a) the phase, 1 — (/3a,s/</ 5, and (b) the magnitude of 
the next-to-nearest tunneling coefficient, 1 — |tiA,s|/|ti|, for 
Xa = 0.001, s = 5 Eh. Results calculated using the MLWFs. 


s - of the underlying continuous Hamiltonian. 

Next, we proceed to test the accuracy of the assump¬ 
tions of the SPS approximation (Sec. IIIB) based on the 
tunneling coefficients calculated from the MLWFs. This 
is done in Fig. where we compare the relative devia¬ 
tions from the average values of the phase, 1 — ipa,b/'-P, 
and of the magnitude of the next-to-nearest tunneling 
coefficient, 1 — |tiyi,B|/|ti|, for xa = 0.001, s = 5. This 
figure demonstrates that the maximum relative deviation 
in both cases is below ^1%. We have verified that this 
holds for all values of s and xa considered here, thus 
justifying the assumptions of the SPS approximation in 
the whole range of parameters. Apart from the relative 
deviation. Fig. |^a) reveals that (pA and (pB exchange 
roles at a ~ 1.5, around the point where the phase gets 
its maximum value, see Fig. id). 


IV. TOPOLOGICAL PHASE DIAGRAM 


(37)-(40l and (41|-(44) exchange roles at the boundaries 
between normal and topological insulator regimes; this 
feature is particularly noticeable in Figs. ia) and id). 
To put in other words, the solution of one set of equa¬ 
tions on one side represents a smooth continuation of the 
solution of the other set of equations in the other side, 
and viceversa. Provided that one chooses the right so¬ 
lution, the calculated values agree very well with those 
of the MLWFs, as already said. In addition. Fig. id) 
reveals an extremely important feature that was absent 
in the original Haldane model: the phase p is limited by 
a maximal value. This behavior, that was already found 
in the parity-symmetric case |20j , implies that p can only 
access a restricted range of values, therefore limiting the 
physically accessible region of the phase diagram. This 
feature will be crucial for the analysis presented in the 
next section, where we shall redraw the topological phase 
diagram in terms of the physical parameters - a, xa and 


The topological state of a system is characterized by 
the so-called Chern number or topological index [31] 

« occ 

C=7r dk'^{dku„k\ X \dku^k ), (45) 

Jbz 


with Uvkif) = e.~^’^"^'ipi,k{r) being the periodic part of 
the Bloch eigenfunctions. Since the band structure of 
the Haldane model consists on a valence and a conduc¬ 
tion band, only the lower energy band enters the sum 
over occupied states in Eq. (45). In order to efficiently 
calculate the Chern number, one can rewrite the expres¬ 
sion in (45) as 



/ dk n{k), 

IBZ 


(46) 


where n(k) stands for the Berry curvature [32]. This 
quantity can be accurately computed by means of the 
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FIG. 5. (Color online) Topological phase diagram of the Hal¬ 
dane model as a function of ip and e/|ti|. The main figure is a 
zoom for tp € [—0.45,0.45], while the inset illustrates the full 
nominal diagram, with p £ [—7r,7r]. Large (green), medium- 
size (red) and small (blue) dots correspond to non-zero Chern 
numbers calculated ab initio for s = 5, 7 and 9, respectively. 
The sign of the Chern number is equal to the sign of the 
phase. The solid (black) line denotes the analytical bound¬ 
ary e/|ti| = 3\/3sin(/9. The vertical dashed lines delimit the 
physically accessible regions. 


Wannier interpolation technique, as discussed in [271 [Ml 
[33] ■ In our calculations, we find that a fine 5000 x 5000 
fe-mesh is required in order to converge the integral of 
Eq. (ra. 


The Chern number represents a topological property 
and takes only integer numbers |31j . Its value is inti¬ 
mately connected to the band structure and the gaps 
opened by symmetry breaking at the Dirac points. If a 
gap is opened solely by inversion symmetry breaking, the 
state of the system is topologically trivial with (7 = 0 . 
On the other hand, if the gap is opened by time-reversal 
symmetry breaking, then the system is found in a topo¬ 
logically non-trivial state with (7 yf 0. When both sym¬ 
metries are broken, the topological state of the system 
depends on the relative strength of the inversion and 
time-reversal symmetry breaking. 

The topological phase diagram of the Haldane model 
has been traditionally drawn as a function of p and 
e/iHi mm- In order to facilitate the discussion, let 


us rewrite here the analytic expression in Eq. (26 1 that 


defines the boundary between the different insulating re¬ 
gions, namely 


1 —j- = ±3-\/3 sin (/?. 
Til 


(47) 


In the original formulation, in which the dependence 
of on a is derived by means of the Peierls substitu¬ 
tion HKH], the whole phase diagram is accessible. How¬ 
ever, since the Peierls substitution is incorrect [20], the 
possible values of p are actually limited to a finite range 
that depends on s, as discussed in section HIC (see e.g. 



FIG. 6. (Color online) Topological phase diagram of the con¬ 
tinuous Hamiltonian in Eq. Q, as a function of a and xa, 
for three different values of the scalar potential amplitude s. 
The non-trivial topological state is indicated by big (green) 
dots for s = 5, medium (red) dots for s = 7 and small (blue) 
dots for s = 9. The black dashed lines represent a guide to 
the eye for the phase boundaries for each value of s. 


Fig. [^d)). This is shown in Fig. where the accessible 
region for each value of s is represented by the vertical 
(dashed) lines. Actually, only a small portion of the nom¬ 
inal phase diagram can be accessed (see the inset), as the 
maximum allowed values of p are much smaller than tt. 
In the figure, the dots represent a non-trivial topological 
state with (7 = ±1. The fact that almost all these points 
lie in between the black solid lines [M] proves that - in the 
allowed accessible region - the phase diagram of the mi¬ 
croscopic Hamiltonian is well described by the analytical 
expression of Eq. (47) for the Haldane model. 


Owing to the above analysis, we suggest that a more 
appropriate way to draw the topological phase diagram is 
in terms of the physical parameters that characterize the 
underlying continuous Hamiltonian, namely a, xa and 
s. This is shown in Fig. [^ where we plot the phase 
diagram in the a — XA plane, for three different values of 
s. Importantly, the Fig. evidences that the topological 
insulating phase with (7 7 ^ 0 shrinks dramatically as the 
system becomes more and more tight-binding (that is, by 
increasing s). Notice that the sign of the Chern number 
in the topological insulator phase {C = ±1) is consistent 
with the sign of a, and independent on the sign of xa- 
Notice also that the probability of finding the system in 
the topological insulator phase increases consistently by 
decreasing the value of |xa|- 

As previously anticipated, the structure of the phase 
diagram is intimately connected to the behavior of the 
gaps at the Dirac points. This is illustrated in Fig. [^ 
where we plot the gaps 5+ and S- as a function of a > 0 
and three different values of xa > 0 for fixed s = 5. Note¬ 
worthy, the gap closing does not take place exactly at 
but in a close-by non-high-symmetry point. The origin of 


this feature has already been discussed in Sec. HC For 








































































FIG. 7. (Color online) behavior of the gaps (squares, solid 
line) and 5_ (triangles, dashed lines) as a function of q, for 
s = 5. The three panels correspond to xa = 2 • 10“"* (a), 
10“^ (b), and 1.9 • 10“^ (c). The latter corresponds to the 
maximal value of \xa\ for which the system can be in the 
topological insulating phase. Grey shaded area corresponds to 
the region where the system is a topological insulator {C = 1), 
whereas the white background identifies a normal insulating 
state ((7 = 0). The inset in panel (b) shows the behavior of 
the gap S- at fc)) (dashed blue line) and at ~ fc)) + 1.68 • 
10“^(l/2, %/3/2)fc_L (red continuous line), around the point 
a « 1 . 8 . 


the specific value s = 5 analyzed here, our calculations 
identify this point at k~ ~ fc)^1.68- 10“^(l/2, y/3/2)kL, 
as shown in the inset of Fig.^b) [35]. We find that the 
gap closing point is slightly shifted for different values of 
s, but lies always very close to k~jj. In all cases, the de¬ 
viation from k'^ represents a minor correction, and can 
be safely ignored in the following discussion. 


Notably, Fig. reveals that the gap has a maximum 
at a ~ 1.0 fci, impliying that the effect of the vector 
potential in opening the gap is limited, as expected from 
Eq. (26). It is also noteworthy that when xa is rela¬ 
tively small, as in Figs. and [^, the gap S- vanishes 
for two different values of a (the role of S+ and S- is 
exchanged for a < 0) . In fact, owing to the non mono¬ 
tonic behavior of as a function of a, see Fig. w and 
Ref. [3D] , there are two different values of a for which Eq. 
(47) can be satisfied (notice that the two values of (/? at 
the phase boundaries may be slightly different due to the 
fact that ti also depends on a). The intermediate region 
between these two values, which is represented by a grey 
shaded area in the figures, corresponds to a topological 
non-trivial state {C = 1) where the effect of time-reversal 
symmetry breaking is stronger than inversion symmetry 
breaking. As mentioned above, the smaller the value of 
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FIG. 8. (Golor online) Gap created at the Dirac points by 
time-reversal symmetry breaking, for three different values of 
the scalar potential amplitude s, fixed value xa = 2 • 10~^. 
Solid (squares) and dashed lines (triangles) denote <5+ and S -, 
respectively. The vertical dashed-dot-dot lines denote the gap 
closing points for the two lowest values of s. 


XA, the larger the region with (7 0 as a function of 

a. By increasing xa, the topological insulating phase 
shrinks and eventually disappears, as shown in Fig. [^c). 

To conclude our analysis, let us discuss why the phase 
diagram of Fig. [^shrinks as s is increased. For such pur¬ 
pose, in Fig. I^we illustrate the evolution of the gap as 
a function of a and s for fixed xa- This figure evidences 
that the maximum of the gap decreases as s is increased; 
in other words, the relative effect of time-reversal sym¬ 
metry breaking decreases with increasing s. As a conse¬ 
quence, even relatively low values of xa can avoid gap 
closing provided s is large enough, as in the case of s = 9 
in Fig. 1^ This, in turn, implies that the phase tran¬ 
sition to the topological insulator phase is restricted to 
smaller values of xa as s is increased, in agreement with 
the phase diagram of Fig. 

V. CONCLUSIONS 

In summary, we have presented an ab initio analysis 
of a continuous Hamiltonian |I8j that maps into the cel¬ 
ebrated Haldane model [T] . The tunnelling coefficients of 
the tight-binding model have been computed by means of 
two independent methods, one based on the maximally 
localized Wannier functions and the other on a closed 
set of analytical expressions in terms of the energy spec¬ 
trum at selected high symmetry points in the BZ. The 
two approaches present a remarkable agreement. In par¬ 
ticular, we have shown that the gaps created either by 
inversion or time-reversal symmetry breaking are very 
well described by the tight-binding model, which repro¬ 
duces accurately the exact behavior. In addition, we have 
calculated the topological phase diagram in terms of the 
physical parameters entering the microscopic Hamilto¬ 
nian, finding that only a small portion of the original 
phase diagram discussed by Haldane can be actually ac¬ 
cessed within this model. Moreover, we have shown that 
the non-trivial topological phase with non-zero Chern 
number is suppressed as the system enters the deep tight- 
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binding regime. We believe that, besides its conceptual 
implications, this work is relevant for a possible exper¬ 
imental implementation of the Haldane model following 
the proposal in Ref. m- 
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Appendix A: Spread of the MLWFs 

Here we analyze the properties of the spread func¬ 
tional of the MLWFs, fl = ES], as 

the amplitude a of the vector potential is varied and the 
system crosses the topological phase boundary. Marzari 
and Vanderbilt showed that this functional can be di¬ 
vided into three parts, namely fl = Qj + flo + floD [23]. 
The term fl/ is gauge-invariant (namely, independent of 
the choice of the unitary transformations (fe) in Eq. 
(®), whereas the diagonal term Qo and the off-diagonal 
term Hod do depend on the gauge choice. In Fig. 
we show the behavior of the three terms of the spread 
as a function of a, for fixed values of s and XA- Here, 
the non-trivial topological phase is indicated by the grey 
shaded area. All the components of the spread show a 
continuous behavior, even across the boundary between 
the trivial and non-trivial topological states. Then, it is 
interesting to note that, while the gauge-invariant term 
rij shows a monotonic decrease as a function of a, the 
gauge-dependent terms flu and floD show a non mono¬ 
tonic behavior that is reminiscent of what we observed 
for the gap (see Fig. and for the complex phase of the 
next-to-nearest tunnelling coefficient in Fig. id). 

We notice that the smooth behavior of the spread 
shown by our calculations differs from an earlier anal¬ 
ysis of MLWFs in the context of the Haldane model per¬ 
formed by Thonhauser and Vanderbilt [36]. There the 
authors found a breakdown of the usual procedures to 
build MLWFs as the system approaches the topological 



FIG. 9. (Color online). Spread of the MLWFs as a function 
of a for fixed values s = 5, xa = 1 • 10“^. The spread is 
decomposed into its gauge-invariant (H/), band diagonal (Ho) 
and band off-diagonal (Hoo) terms. Note the 10® factor in 
the case of Ho and Hoo- 

phase boundary, resulting in a divergence of the spread 
functional. The fundamental difference between our ap¬ 
proach and the one followed in Ref. El] resides in the 
set of bands considered for the construction of the ML¬ 
WFs. In fact, whereas our set includes both the valence 
and conduction bands, their approach included only the 
valence band. This is a crucial difference, since the net 
Chern number of a single band in the topological phase is 
finite, therefore it becomes impossible to choose a smooth 
periodic fe-space gauge of the Bloch orbitals and the pro¬ 
cedure for constructing the MLWFs fails. In our case, 
in contrast, the net sum of the Ghern numbers of the 
valence and conduction bands remains null, hence there 
is no formal impediment for the construction of the ML¬ 
WFs. 


Appendix B: Numerical calculation of the spectrum 

Both the calculation of the exact Bloch spectrum of the 
continuous Hamiltonian of Eq. Q and the construction 
of the MLWFs require a standard Fourier decomposition 
that here is adapted to account for the presence of the 
vector potential. We express the eigenstates ipnkir) of 
the Hamiltonian as 

tpnk{r) = ^ (Bl) 

G 

with G the reciprocal vectors and Cnk+G the expansion 
coefficients. The vector potential acts as A{r) ■ p + p ■ 
A{r) = —2ihA{r) ■ V^, introducing a non-local term 
when acting upon an eigenstate iprikif): 

iA{r) ■ Vr^nkir) = -A{r) ■ ^ Gcnk+Ge^^ ^ (B2) 

G 

Numerically, we found that a large number of G vectors 
are needed in order to converge the above term due to the 
presence of the gradient. In particular, the above term 
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requires an energy cutoff of 50 En, whereas the rest of 
the terms in the Hamiltonian are converged with 10 En. 

Finally, for extracting the tight-binding parameters us¬ 
ing the formulas discussed in section [Til B[ we have used 
a direct diagonalization of iJo in Eq. ([^ by means of a 
standard Fourier decomp ositi on. In this case, the vec¬ 
tor potential term in Eq. (B2) is transformed into a non 


diagonal matrix in momentum space. 
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